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Abstract 

We consider the problem of minimizing a function represented as a sum of submodular 
terms. We assume each term allows an efficient computation of exchange capacities. This holds, 
for example, for terms depending on a small number of variables, or for certain cardinality- 
dependent terms. 

■ A naive application of submodular minimization algorithms would not exploit the existence 

^"*^ | of specialized exchange capacity subroutines for individual terms. To overcome this, we cast 

^ ■ the problem as a submodular flow (SF) problem in an auxiliary graph, and show that applying 

U | most existing SF algorithms would rely only on these subroutines. 

We then explore in more detail Iwata's capacity scaling approach for submodular flows [T5] . 
In particular, we show how to improve its complexity in the case when the function contains 
• cardinality-dependent terms. 

o ■ 

On 

1 Introduction 

In this paper we consider the problem of minimizing an objective function of the following form: 
O ' 

f(S)= E fQ( Sr] Q) V5 ^ V (!) 

Qe<2° 



X 



Here V is a set of nodes, Q° C 2 V is a set of subsets of V, and /q : 2^ — > R are submodular 
functions. 

Function / is itself submodular, and thus can be minimized in polynomial time. The current 
fastest strongly polynomial algorithms are those of Orlin [23] and Iwata-Orlin [21], which take 
time 0(n 5 EO + n 6 ), where n = \V\ and EO is the time to run the value oracle for f(S). The 
fastest weakly polynomial algorithms are those of Iwata [19] and Iwata-Orlin |21] which run in 
time 6(n 4 EO + n 5 ). 

However, applying a general-purpose submodular minimization algorithm may not be the most 
efficient technique, since it does not exploit the special structure of /. It is often the case that 
terms /q have a special form that allow an efficient computation of exchange capacities, which 
are defined in the next section. Roughly speaking, this means that we can efficiently minimize 
function /q(S) — z(S) for any vector z € MP . (As usual, z(S) denotes ^ ig g Zi.) The main goal 
of this paper is to develop an algorithm that can exploit the existence of specialized exchange 
capacities subroutines. 

To achieve this goal, we use the framework of submodular flows (SF) introduced by Edmonds 
and Giles [8]. We show that the problem of minimizing / can be cast as a particular SF instance 
in an auxiliary graph, so that computing exchange capacities for the new problem is equivalent to 
computing exchange capacities for individual terms /q. Most existing algorithms for submodular 
flows rely on the exchange capacity oracle, which gives the desired result. 
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We then present a capacity scaling technique for solving the problem. Its complexity is 0((n + 
a o){ n + Pq) l°g U) where U is an upper bound on function values and ckq, /3q depend on 
the type of term /q : 

(a) If |Q| = 2then (olq, /3q) = (1, 1). 

(b) If /q(5) = ff(|5|) then (a Q ,f3 Q ) = (\Q\, \Q\). Note, g(-) must be concave. 

(c) If f Q (S) = g(\S n Q'\, \S n Q"|) where Q',Q" are disjoint subsets of Q then (a Q ,/3 Q ) = 
(IQIMQD- 

(d) For any other term /q we have (cxq,/3q) = (|Q| 2 ,|Q| 2 + |Q| • ho) where Hq is the time of the 
exchange capacity oracle for the (scaled version of) /q. 

In (b) and (c) we assume that function g can be evaluated in O(l) time. For cases (c) and (d) we 
use the scaling technique of Iwata [18] . 

Applications Functions with terms of the form (a)-(c) have recently appeared in computer 
vision applications. Terms (a) and (b) were used for the image segmentation problem (22|, 125] . 
while terms (a) and (c) were used for co-segmenting two images containing a similar object |17| . 
(The latter work used terms of the form Jq(S) = — c ■ \S fl Q'\ ■ \S D Q"\ with c > 0.) 

Note, objective functions used in computer vision very often have form (pQ) where \Q\ is quite 
small (2,3,. . .). Terms /q encode interactions between neighboring pixels. Currently, researchers 
restrict themselves to functions that can be reduced to a minimum s-t cut problem (see discussion 
in |26j). since minimizing general submodular functions is too expensive in practice. Our work 
may remove such restriction. 

Related work The problem of minimizing functions of the form (pQ) was studied by Cooper [5], 
who formulated a linear program with an exponential number of constraints and showed that its 
optimal value coincides with the minimum of /. The formulation that we will use closely resembles 
that in [5]. Note, however, that the question of how to solve this LP efficiently was not addressed 
in [5], and a connection to the submodular flow problem was not given. 

It is known that in certain cases the problem can be reduced to a minimum s-t cut problem in 
a graph with auxiliary nodes. Billionnet and Minoux [2] showed that this can be done for functions 
with cubic terms, i.e. when \Q\ < 3 for all terms /q. Reductions for certain subclasses with higher 
order terms were given by Freedman and Drineas |10| . Kohli et al. |22j and Zivny and Jeavons |27| . 
The resulting maxflow problem could be solved e.g. in 0(min(n 2 / 3 , m 1 / 2 )mlog(h 2 /m) log U) time 
by the algorithm of Goldberg and Rao [16], where n, rh are the number of nodes and edges in the 
constructed graph and U is a bound on edge capacities. 

On the negative side, Zivny et al. [26] proved that some submodular terms with \Q\ = 4 do 
not admit such a reduction. Even if the reduction exists, it may result in a graph which would 
be prohibitively large in practice. Consider, for example, terms of the form /q(S) = g(\S\) where 
g is concave. The reduction of Kohli et al. |22] adds b extra nodes and b\Q\ extra edges for each 
term fn, where b is the number of breakpoints of the piecewise-linear concave function g. If g 
is strictly concave (as in the application of [25]) then b = \Q\ — 1, so there would be 0(|Q| 2 ) 
edges. In contrast, our technique uses only 0(|Q|) memory. The same holds for the function 
/ Q (5) = -c-\Sn Q'\ ■ \S n Q"\ used in [32]. 

Fujishige and Iwata |12] considered functions of the form f(S) + g(\S\) on a distributive lattice 
where / is submodular and g is concave. It was shown that the problem is equivalent to a 
parametric problem: minimize function of form f(S) + c\(S) for all values of A, where {c\}\ is a 
certain family of non-increasing vectors in MX . 
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2 Problem formulation 

Let Q be the set obtained from Q° by removing all singleton subsets of the form {i}, i E V. Thus, 
\Q\ > 2 for all Q E Q. Without loss of generality we assume that function / is given by 

f(S) = Y,dt+ £ c S4 +^/ Q (5nQ) (2) 
ies iev-s QeQ 

where ca, c s { are non- negative numbers and each term /q satisfies the following condition^] 

mmf Q (S) = f Q (0) = f Q (Q) = O (3) 

SC2 ( * 

Base polyhedron and exchange capacities The base polyhedron [7] of /q is defined as 

B(f Q ) = W Q £R Q \ VQ (S)<f Q (S) V5CQ, <p Q (Q) = f Q (Q) = 0} (4) 

Given a vector </?q E B(/q) and distinct nodes i, j € the exchange capacity CQij is the maximum 
value of e > such that the operation ipQi := tpQi + e, (fQj := ifQj — e keeps ipQ in Clearly, 

c Qii = min{/ Q (5) | i E 5 C Q - {j}} , / Q (5) = f Q (S) - VQ (S) (5) 

Computing CQij is equivalent to minimizing a submodular function. This can be done in polynomial 
time by a number of general-purpose submodular minimization algorithms. Furthermore, for many 
choices of /q there exist more efficient specialized techniques. 

A remark on notation: in this paper we always use "bar" (cQij, /q, • • •) to indicate "residual" 
values, i.e. values that take into account current flow. 

Maximum flow formulation Let us construct a directed capacitated graph G = (N, A, c) as 
follows. The set of nodes will be N = {s,t} U V U<g g g Q* where s,t are the source and the sink 
and Q* = {Qi \ i E Q} is a unique copy of Q. Here Qi is a shorthand notation for the pair (Q, i). 
The set of arcs will be 

A = {(i,Qi),(Qi,i)\ieV,QieN}\J{(s,i),(i,s),(i,t),(t,i)\ieV} 

Arc capacities c s j, cu are the same as in ([2]). Arcs to the source and from the sink have zero 
capacity (cj s = cu = 0), and all "internal" arcs have infinite capacity (c^Qj = cq^ = +oo). 

A flow ip is a vector in M" 4 . For a subset Q E Q we denote 9?q E MP to be the vector with 
components ipQi = <p% t Q%- We also denote value{(p) = YieV *° ^ e ^ ne amoun t of flow sent from 
the source. We will consider the following maximum flow problem: 

max value(tp) s.t. (6a) 

Puv = —^Pvu V(it, i>) E A (antisymmetry) (6b) 

tp a < c a Vo G i (capacity constraints) (6c) 

= Vi E V (flow conservation for V) (6d) 

(u,i)eA 

(fQ E B(fo) \/Q E Q (base polyhedron constraints) (6e) 

Note, if ^ is feasible then we also have value(ip) =^ ieV (pa since Yiev fsi-Yiev Pn = J2ievi l Psi + 



1 If term /q with fg(0) = does not satisfy Q then we can replace it with the sum ipq(S) + f(S) where 
f(S) — f(S) — ipci(S) and ipQ is a vector in the base polyhedron of /q, which can easy be computed by a greedy 
algorithm of Edmonds [7] . 



The linear program ([6]) is very similar to that in [5], with some minor differences; for example, 
the "balance" constraint ipq{Q) = is not present in [5]. 

The rest of the paper is organized as follows. Section [3] gives a reduction of problem ([6]) 
to a submodular flow problem, which leads to a number of algorithms for solving ([6]). Section [4] 
describes a pseudo-polynomial augmenting path algorithm, which is a specialization of the standard 
augmenting path algorithm for submodular flows. By analyzing the algorithm we will prove that 
the maximum of © coincides with the minimum of /. Section [5] presents a scaling version of the 
augmenting path algorithm, while section [6] discusses some implementational issues and states the 
complexity of the algorithm. 

The reader may choose to skip the next section; familiarity with the submodular flow problem 
will not be necessary for understanding the augmenting path algorithm. 



3 Reduction to a submodular flow problem 

We will consider a directed capacitated graph G' = (N,A',c) where A' = A U {(s,t), (t, s)} and 
the capacities of the new arcs are q s = +oo, c s t = 0. If ip G R" 4 is a flow in G' and u is a node in 
N then dip{u) = Yl( v u)eA' Pvu will denote the amount of flow that comes into u. 

Let us recall a definition of a submodular flow problem for a graph G' [81 [15] . Assume that 
each arc a G A' has a cost d a , and let g : 2 N — >• R be a submodular function with g(0) = g(N) = 0. 
Then the problem is defined as 

max y^ 3CA , s - t - ( 7a ) 

<Puv = ~<Pvu V(n, v) G A' (7b) 
fa <c a Va G A' (7c) 
dip G B(g) (7d) 

where B(g) is the base polyhedron of g: 

B{g) = {z £ R N | z(X) < g{X) VX C N, z(N) = 0} (8) 

In order to simulate problem Q, we set arc costs as follows: dt s = 1 and d a = for all other 
arcs a. Function g is defined by 

g(X) = fQ( xQ ) 
QeQ 

where we introduced notation X® = {i € Q | Qi G X}. 
Proposition 1. Problems ([6]) and (JT]) are equivalent. 

Proof. Suppose that ip € R" 4 is a feasible flow for problem ([6]). Let us extend it to a flow in G' 
by setting tp ts = value(p), p s t = —value(p). Clearly, conditions (|7bj) and (f7c|) are satisfied. It is 
also easy to check that z = dp E B(g). Indeed, we have Z{ = for i G V U and zq^ = 99Qj 

for Qi € AT. Conditions <^Qj € B(/q) then imply that z(iV) = and for any X Q N there holds 
z(X) = EQes^(^ Q ) < EQg S /Q(^ Q ) = ffPO- Thus > V is a feasible flow for problem ©. 
Furthermore, the values of objective functions of ([6]) and (|7|) coincide. 

Conversely, suppose that ip € M. A is a feasible flow for problem ([7]); let us show that its 
restriction to A is feasible for problem ©. Conditions (|6b|) and ([6c| follow from ()7b|) and ([7c| . 
Denote z = If A is a subset of A with g(X) = g(N — X) = then z £ B(g) implies 

z{X) < g(X) = and -z(X) = z(N - X) < g(N - X) = 0, so z(X) = 0. Applying this fact for 
subset X = {i} yields (|6dj) . and applying this fact for subset X = Q* yields constraint pq(Q) = 0, 
which is a part of (f6e|) . Finally, if S C Q then pQi(S) = z(S*) < g(S*) = Jq(S) where we denoted 
S* = {Qi | i G 5}. Thus, ^ Q G 5(/ Q ). □ 
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Exchange capacities Most submodular flow algorithms rely on the following operation: given 
a feasible flow ip G IR" 4 with z = dp G B(g) and distinct nodes u,v G N, compute the exchange 
capacity c uv = mmx{g(X) \ u G X C N — {v}} where g(X) = g{X) — z(X). The proposition 
below shows that computing these capacities is equivalent to computing exchange capacities CQij 
for individual terms /q with respect to flow p (given by eq. ©). 

Proposition 2. c uv = CQij if (u,v) = (Qi,Qj) and c uv = otherwise. 

Proof. As shown above, z% = for i G V U {s, t}, therefore z(X) = X^QeQ PQii-X ) for all subsets 
X Q N. This implies that 

g(X) = ^ / Q (X«) (9) 
QeQ 

The fact that v?q G #(/q) also implies min 5 cQ/Q(5') = /q(0) = /q(Q) = for all Q G Q. 
Therefore, if (u, v) = (Qi, Qj) then the minimization problem minx{^(^) j u G X C — {f }} has 
a minimizer X C Q* } and thus c ra = mmx{g(X) \ u £ X Q Q* — {v}} = c~Qij- Now suppose that 
(u,v) / (Qi,Qj). Let U C N be the "completion" of u: U = {u} if u G V U {s,t} and U = Q* 
if it = Qi. There holds v ^ U since we assumed that (it,t>) ^ (Qi,Qj) and it, w are distinct. We 
have 5(C/) = 0, and thus c uv = 0. □ 

Problem (|7|) is actually a maximum submodular flow problem, which is a special case of the 
more general minimum cost submodular flow problem (see survey [15] ) . The former problem can 
be solved in time 0(\N\ 3 h) by a push-relabel method of Fujishige and Zhang |13j . where h is the 
time of the exchange capacity oracle (see also |20| . section 3.1). Clearly, for certain functions / this 
complexity can be better than bounds 0(n 5 EO + n e ) and 0(n 4 EO + n 5 ) for submodular function 
minimization. 

In our case h is the maximum time of oracles over individual terms. This appears to be a rather 
crude way of estimating the complexity, as it does not take into account the structure of individual 
terms. We conjecture that a more careful analysis of the algorithm can give a bound which better 
illustrates contributions of individual terms. In the subsequent sections we will give an example 
of such a bound for a capacity scaling augmenting path algorithm applied to problem ([6]). 



4 Augmenting path algorithm 

A shortest augmenting path algorithm for a problem equivalent to maximum submodular flows 
was given by Fujishige [H]. We now describe its application to problem ([6]), and prove that the 
value of the maximum flow coincides with the minimum of /. We will generalize the problem 
slightly: we assume that capacities Ci s and cu are non-negative numbers which are not necessarily 
zero. (We will need this in the next section.) 

Given a flow p, the residual capacity for arc a G A is defined as c a = c a — <p a . Similarly, we 
define "residual functions" /q by /q(5) = /q(S) — <Pq(S) for S C Q. It can be seen that if p 
satisfies antisymmetry and conservation constraints (I6bp . (I6dp then for any S (IV there holds 



f(S) = value{p) + Y J Cit+ £ c sl + ^ f Q (S n Q) (10) 



ies i&v-s 



Indeed, subtracting © from (HDD gives Y. i& v <Psi ~ Ei G s fit ~ Hi&vs <Psi ~ EqgS EiesnQ <Pi,Qi = 
^2i£S + + Sqigtv ^QM = 0- All residual values for a feasible p> are non-negative, so 
equation (fTUj) implies the weak duality relationship: 

max{value(p) \ p is feasible } < mm{f(S) \ S C V} (11) 



5 



Given a feasible flow ip, let A be the following set of arcs: 

A = {a e A | c a > 0} [j A Q , A Q = {{Qi, Qj) \ i,j £Q,i^ j, c Qij > 0} (12) 

Proposition 3. If there is no path from s to t in (N, A) then the set S = {i £ V \i is reachable from 
s in (N, A)} satisfies f(S) = value(p>), and therefore (p is a maximum flow and S is a minimizer 
off- 

Proof. It suffices to show that every term in the RHS of (|10p (except maybe for the first term 
value(p)) is zero. If i & S then cu = 0, otherwise t would be reachable from s. If i £ V — S 
then c s i = 0, otherwise i would belong to S. Consider the term for subset Q € Q, and denote 
5' = S (~l Q. For each pair of nodes i € S' , j G Q — S' function /q must have a minimizer Sij 
with i € Sij C Q — {j}, otherwise we would have CQij > so node j could be reached from i via 
arcs (i, Qi), (Qi, Qj), (Qj,j) £ A and thus j would be in S. The submodularity of /q implies that 
the set \J ie s' HjeQ-s' Sij 1S a minimizer of /q as well. The latter set coincides with S' = S n Q, 
therefore f Q {S HQ) = 0. □ 

Now suppose that there exists a path P from s to t; such a path is called an augmenting path. 
Clearly, we can send some flow 5 > along the path! so that the flow would remain feasible and 
value((f) would increase by 5. This leads to 

Proposition 4 (Strong duality). The value of the maximum flow in ([6]) coincides with the mini- 
mum of f. 

Proof. Let (p be a maximum flow. There can be no augmenting path for ip, otherwise p would not 
be maximal. The claim now follows from proposition [3J □ 

From now on, we assume that all capacities c^, c it and values /q(S') for S C Q are integers 
bounded by constant U. A maximum flow can then be computed in pseudo-polynomial time by 
the following augmenting path algorithm: 

50 Set p a = for all arcs a. 

51 Construct set of arcs A as in (|12D . 

52 Find a shortest path P from s to t in (iV, A); if no such P exists, terminate. 

53 Send 1 unit of flow along P and go to step 1. 

Note, it is well-known that for integer-valued submodular flow problems sending 1 unit of flow 
along a shortest augmenting path preserves flow feasibility p3]. In our case we can relax slightly 
the requirement that P is shortest; we only need P to be minimal: 

Definition 5. Let P be a simple (i.e. node- disjoint) path in (N,A). We call P minimal (with 
respect to (N,A)) if the following property holds: if (Qi,Qj), (Qi',Qj') are two distinct arcs in 
the path (occurring in this order) then A does not have arc (Qi,Qj'). 

Clearly, any shortest augmenting path from s to t is minimal. In Appendix A we prove that 
sending one unit of flow from s to t along a minimal path preserves flow feasibility. 

It is not difficult to show that sets Aq are transitive, i.e. (i,j), (j, k) S Aq implies (i, k) € Aq 
(see Appendix A). Thus, if P is minimal then (Qi, Qj) € P implies that the previous arc in P is 
(i, Qi) and the next arc is (Qj,j). The operation of sending flow through these three arcs will be 
referred to as "sending flow from i to j via Q" . 

2 Sending flow S along arc (u, v) £ A denotes the operation <p uv := <p uv + S, ip vu := ip vu — 8. Sending flow 8 along 
arc {Qi, Qj) £ Aq does not change ip. 



G 



so 



For each Q £ Q set ipQ := AdjustFloWg^Q) to make sure that ipQ £ B(/q). Ad- 
just other flow components so that 99 satisfies antisymmetry and flow conservation 
constraints: 



• Set ip Qit i := -ipi,Qi for all Qi £ N. 

• For each node i £ V compute 5 = Yli(ui)&A fui] if <5 > 0, send 5 units of flow 
back to the source via arc (i, s), otherwise send —5 units of flow from the sink 
via arc (t, i). 



SI 



Construct set of arcs A 



as follows: 



{{u,v) £ A 



c U v > r A i } U 



(13) 



52 Find minimal path P in (iV, A ); if no such path exists, terminate. 

53 Send [A] units of flow along P and go to step SI. 

Figure 1: A-phase. Definitions of function fg, procedure AdjustFlowg(c/?) and set A^ for dif- 
ferent types of terms /q are given in sections I5.1H5.31 

5 Capacity scaling algorithm 

We now apply a scaling technique to get a weakly-polynomial algorithm. As usual, the algorithm 
works in phases. Each phase is associated with a number A = 2 l , I = —1, 0, 1, 2, . . .; we call it a 
A-phase. To initialize, we set A = 2r i °S2 c/ l and cp a = for all arcs a £ A. After completing the 
A-phase we divide A by 2 and proceed to the next phase (or terminate, if A = 1/2). The A-phase 
is described in Figure [1] This description uses the following yet undefined objects: 

• fg is a submodular function. When A = ^, function fg coincides with /q. 

• AdjustFloWg (<pq) is a procedure that outputs a vector in B(fg) whose components are 
integer multiples of [A] . 

• Aq is a subset of arcs of the form (Qi,Qj) where i,j are distinct nodes in Q. Set Aq is 
transitive, i.e. (Qi,Qj),{Qj,Qk) £ Aq for distinct k £ Q implies (Qi,Qk) £ Aq. When 
A = |, set Aq coincides with the set Aq defined in (fl"2"l) . 

Definitions of these three objects will depend on the type of term /q; different cases are considered 
in sections 15.1115.31 Set Aq will be defined in such a way that each augmentation keeps flow ipQ in 



It is clear that each A-phase maintains the following invariants: (i) components of flow tp are 
integer multiples of [A]; (ii) cp is a feasible A-flow, i.e. it satisfies antisymmetry (I6b|) . capacity 
(f6c|) . flow conservation constraints (|6dp . as well as base polyhedron constraints cpQ £ B(fg). (We 
assume that capacities Cj s , cu are infinite, so that sending flow to the source or from the sink in 
step SO is always feasible). To estimate the complexity, we will use values q.q (to be defined in 
sections l5.1H5.3|) that satisfy 



where (p° is the flow in the beginning of A-phase, S is the set of nodes in Q reachable from s 
in the graph (N, A 2 ^) constructed with respect to flow (p°, ip = AdjustFlowg((/9°) and fn(S) = 





(14) 



ieQ 
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fn(S) — ipq(S). Values ag can be used for estimating the number of augmentations (a proof is 
given in Appendix B): 

Proposition 6. Each A-phase terminates after at most 2n + X^Qeg a Q augmentations, and so 
the whole algorithm performs 0((2n + ^2q £ qCyq) log U) augmentations. 

To complete the description of the algorithm, we need to provide constructions for different 
types of terms /q. In sections I5.1II5.3I below we consider three types: pairwise terms, cardinality- 
dependent terms and general terms. 

5.1 Pairwise terms 

First, we consider the case when \Q\ = 2, which occurs very frequently in applications (see e.g. [3] 
for a survey of applications in computer vision). We define f§ = f Q for all A. This means that 
procedure AdjustFlowg(y?Q) can simply return ifQ - it is guaranteed to belong to B(fq) = S(/q A ). 
Let Q = Constraint ipQ G B(fq) can be written as 

<PQi < /$({*}) tpQj < /Q({j}) <PQi + <PQj = (15) 

The set of arcs Aq is constructed as follows: we add arc if CQij = /q({^}) > an d arc 
(j, i) if CQji = /g({i}) > [A]. Clearly, we can always push [A] units of flow through the added 
arcs - constraints (|15p will be preserved. 

It is easy to see that we can take uq = 2. Indeed, let S be the set used in eq. (fl4j) . Since 
(p = (p°, we need to show that f Q (S) < 2 [A]. If S = {i} then (Qi,Qj) (£ A 2A , therefore 
Iq(S) < \2A] - 1 < 2 [A] . The case S = {j} is similar. If S is empty or equals Q then f Q (S) = 0. 

5.2 Cardinality-dependent terms 

Let us now assume that /q(5*) for S C Q depends only on \S\. Thus, /q(5) = 9(15*1) where 
g is a concave function. As above, we define fg = fQ for all A, and accordingly procedure 
AdjustFlow((^Q) simply returns ifQ. Below we describe how to construct set Aq. 

For integer numbers a < b let [a..b] be the set of integers in [a, b]. We can assume that g(-) 
is defined only on [0..m] where m = \Q\. We denote z = ifQ, so zi = fi t Qi = —y>Qi,i for i € Q. 
For a vector z € MP we also denote (z 1 , . . . , z m ) to be the sequence of values of Z{ sorted in the 
non-increasing order. Thus, z k is the k-th largest number among values z i: i G Q. For a node 
i € Q define 

L(i) = mm{k £ [l..m] \ z k = Zi} , R(i) = max{/c G [l..m] \ z k = Zi] 

Let us define "residual" function g(-) by g(k) = min{/Q(5) | 5 C Q, \S\ = k} for k € [0..m]. 
We have k 

g(k) = g(k) - max{z(5) | S C Q, \S\ = k} = g(k) - £ z fc ' (16) 

fc'=i 

Clearly, constraint z G B{fo) is equivalent to the following conditions: (i) function g{-) is non- 
negative, i.e. g(k) > for all k G [0..m]; (ii) z(Q) = 0. 

Recall that sending flow [A] from i to j via Q denotes the following operation: Zi := Zi + 
[A] , Zj := Zj — [A] . Next, we describe the effect of this operation on function g(-). Three cases are 
possible (we assume that we are in a A-phase, so all components of vector z are integer multiples 
of \A]): 

• Zi < Zj — 2 [A] . The change in the sequence (z 1 , . . . , z m ) is 

(...,o,-rA],o,...,o,+rA],o,...) 

where — [A] is in the position R{j) and +[A] is in the position L(i). Therefore, the effect 
of the operation is that all values g(k) for k G [R(j)..L(i) — 1] are increased by [A]. 
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• Zi = Zj — [A]. The values Zi and Zj are swapped, therefore the sequence (z 1 ,. 
function g(-) do not change. 

• Zi > Zj. The change in the sequence (z l , . . . , z m ) is 

(...,0,+[Al,0,...,0,-rAl,0,...) 

where + |~A] is in the position L(i) and — [A] is in the position R{j). Therefore, all values 
g{k) for k € [L(i)..R(j) — 1] are decreased by [A]. 

In the first two cases function g(-) cannot become negative, thus sending [A] units of flow from 
i to j via Q is always possible if Z{ < zj. Accordingly, we add arcs (i, j) to Aq for all pairs of nodes 
i,j € Q with Zi < zj. If Zi > Zj then we can send flow if and only if miiike[L(i)..R(j)-i] 9~(k) > [A]. 
However, if we add all arcs that satisfy this constraint then sending [A] units of flow through 
multiple arcs of Q along a minimal path could make some values g(k) negative. To prevent this, 
we add to Aq those arcs with Zi > Zj that satisfy the following constraint: 

min g(k) > 3A/2 (17) 
ke[L(i)..R(j)-i] 

The proposition below shows the correctness of this construction, and gives a bound on an. A 
proof is given in Appendix C. 

Proposition 7. (a) Set Aq is transitive, (b) Sending [A] units of flow through a minimal path 
P in (N,A A ) preserves constraint ipg = z £ B(fn). (c) Eq. (Tn|l is satisfied by qq = 3(m — 1). 



5.3 General submodular terms 

For general terms we can use the technique of Iwata |18j . /q is defined as 

f§(S) = A • [f Q (S)/A\ + [AJ • b(S) VS C Q (18) 

where 6(5) = |5| • \Q — S\. As shown in |18j . this function is submodular. The set Aq includes all 
arcs (Qi,Qj) that have non-zero residual capacity with respect to function fg(S) = fg — (pn(S). 
Clearly, values of fg(S) are integer multiples of [A], so results in section H] imply that pushing 
[A] of flow through a minimal path in (N,A A ) preserves constraint ifQ € B(fg). 

Procedure AdjustFlowg^g) works as follows. First, define vector tp'g by ip'g i = — m[A] 
where m = \Q\. Vector ip'g belongs to submodular polyhedron 

P(f£) = {<PQeR Q \<P Q (S)<f£(S) VSCQ} (19) 

Indeed, for any S C Q we have ipf Q (S) = f° Q (S) - m\A] ■ \S\ < f^ A (S) - m[A] • \S\ < f§(S). 
Since (p'g G P(fg), there exists vector ipg € B(fg) with ip'g > ipg, which can be found by 
a greedy algorithm starting from (p'g Theorem 3.19]. This ifQ is taken as the output of 
AdjustFlow^c^g). 

It can be seen that cxq = 0(m 2 ). This follows from three facts: (1) fg A (S) = where S is the 
set used in eq. (jTlj) and fg A is the residual function with respect to flow ip°; (2) \ fg \S) — fg(S)\ = 
0(m 2 [A]); (3) ZieQ 1^ < ™> [A] +J2 ieQ \<PQi -v' Qi \ = ™ 2 [A] + (p Q (Q) -cp'g(Q) = 2m 2 [A] . 

Note, procedure AdjustFlowg(<^Q) used in [18] is slightly more complicated; in particular it 
takes into account set S used in (|14p . However, both techniques lead to ug = 0(m 2 ). 
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6 Efficient implementation 



We now discuss how implement steps SI and S2 of the algorithm, i.e. how to find a minimal 
augmenting path. Set A A contains 0(n + \Q\ 2 ) arcs, so a naive computation would take 
0(n + \Q\ 2 ) time. However, this can easily be improved: it can be seen that an explicit 
construction of A^ is not required. 

We will use a breadth-first search (BFS) for computing a shortest path from s to t in (JV, A^). 
Each node Qi G N will have flag REACHED((5«), which is set to false at the beginning of BFS. 
We assume that each term /q supports operation GetNeighborsg^i) for a node Qi G N with 
REACHED(<5i) = false. This operation is defined as follows: 

• Compute S = {Qj G N \ (Qi, Qj) G Aq, REACHED (Qj) = false}. 

• Set REACHED(Qj) := true for Q j G S U {Qi}. 

• Return £ as a linked list. 

Flags REACHED(Qi) will not be modified by any other operation (except that they are reset to 
false at the beginning of BFS). 

It is straightforward to implement the BFS procedure using operations 
GetNeighborsg(<5i). The running time of one augmentation (steps S1-S3) will then be 0(n + 
SqgS Pq) where (3q for a fixed Q G Q is the combined time taken by calls to GetMeighborsg(Qi), 
plus the time for sending flow through Q in step S3 (which may update internal structures for 
Q). In Appendix D we show how to implement GetNeighborsq^i) so that /3q = 0(|Q|) in the 
following cases: 

• f Q (S) =g(\S\) for SCQ. 

• f Q (S) = g(\S n Q'\, \S n Q"\) where Q',Q" are disjoint subsets of Q. 

The second case relies on the algorithm of Aggarwal et al. pQ which computes row minima of a 
totally monotone matrix in linear time. For a general submodular term /q a naive implementation 
of GetNeighborsg(Qi) would make \Q\ — 1 calls to the exchange capacity oracle for fg, giving 
(3q = 0(\Q\ 2 h,Q) where \iq is the oracle's complexity. However, the set {(Qi,Qj) \ (Qi,Qj) G Aq} 
can alternatively be obtained from the minimal minimizer in argmin{/g (S) \ i G S C Q}. It 
is natural to assume that computing such minimal minimizer also takes time Hq. Under this 
assumption (3q = 0(\Q\ 2 + \Q\ ■ Hq). Combined with proposition [6j this leads to the overall 
complexity stated in the introduction. 

7 Conclusions and future work 

In recent years there has been an increased interest in the computer vision community in using 
submodular functions of the form ([T]) with high-order terms [22l [23J EH [171 [6] . So far, researchers 
restricted themselves to functionals that can be transformed to pairwise terms by introducing 
auxiliary variables. The main goal of this paper is to advocate a more direct approach which could 
extend the set of practically tractable functionals. 

To our knowledge, our bound 0({n+^2,Q oiQ)(n+^2,Q Pq) log U) is the first one for minimization 
problem ([T]) that shows contributions of individual terms. It is quite likely, however, that it can 
be improved further. Indeed, the capacity scaling algorithm of Iwata [18] that we built on is not 
a state-of-the-art. In the future we plan to investigate applications of alternative submodular flow 
algorithms, such as the capacity scaling algorithm of Fleischer et al. [9] that improves on [18], or 
the push-relabel method of Fujishige and Zhang |13j . 



10 



Appendix A: Minimal augmenting paths 

First, let us show the set Aq defined in fjl2[> is transitive, i.e. if i,j, k are distinct nodes in Q then 
(Qi,Qj), (Qj,Qk) G Aq implies (Qi,Qk) G Aq. Suppose not, then CQik = 0. This means that 
/q(S) = for some subset S with i G S C Q — If i G S 1 then CQj^ = and (Qj, Qk) £ Aq, 
and if j £ S then CQij = and (Qi, Qj) ^ - a contradiction. 

Assume that the problem is integer- valued. It is straightforward to check that sending one 
unit of flow along a minimal path in (JV, A) from s to t preserves antisymmetry (|6b|) . capacity ([6c]) 
and flow conservation (|6d|) constraints. We now prove that if P is a minimal path in (N, A) 
whose endpoints belong to V then sending one unit of flow along P preserves base polyhedron 
constraints (|6el) . Note, P is not an augmenting path: it does not go from s to t. However, the 
operation of sending flow along P and the minimality of P are still well-defined. 

We use induction on the length of P. If P is empty then the claim is trivial. Suppose P is 
not empty; since P is minimal and Aq are transitive, P must have the form P = P1P2 where 
Pi = ((i, Qi), (Qi, Qj), (Qj,j)) and i,j are distinct nodes in Q G Q. Since (Qi, Qj) G Aq, sending 
one unit of flow along P\ preserves base polyhedron constraints. We prove below that after 
sending this flow P2 remains a minimal path in (JV, A); the claim will then follow by the induction 
hypothesis. 

Clearly, we need to consider only arcs in Aq - subsets Aq> for Q' G Q — {Q} are not affected. 
Let us denote /q to be the residual function after sending the flow and Aq to be the corresponding 
set of arcs. We have /q(S) = /q(S) — [i G S] + [j G" S] for 5 C Q, where [•] is the Iverson bracket: 
it is 1 if its argument is true, and otherwise. We need to show two facts: 

(a) if (Qk,Ql) G P2 then (Qk,Ql) remains in Aq; 

(b) if (Qk, QI) and (Qk', QI') are two distinct arcs in P2 occuring in this order then arc (Qk, QI') 
still does not belong to Aq. 

Proof of claim (a) For a set S C Q denote [S] = ([i G S],[j G S], [k G S],[l G S}). 
If the claim is false then there exists S with [S] = (?, ?, 1, 0) and Jq(S) = 0. Since (Qk,Ql) G Aq be- 
fore sending the flow, we must have Jq(S) = }q(S) + [i G S] - [j S] > 0, therefore [S] = (1, 0, 1, 0) 
and Jq(S) = 1. By minimality of P arc (Qi, QI) was not in Aq before sending the flow, therefore 
there exists another set S' with [S'] = (1,?,?,0) and Jq(S') = 0. Since (Qi,Qj),(Qk,Ql) G Aq 
we must have [S'\ = (1, 1,0,0). 

By submodularity f Q (SnS') + f Q (SuS') < f Q (S) + f Q (S') = 1, so one of the set SnS', SUS' 
is a minimizer of /. We have [SnS'] = (1,0,0,0) and [5 US'] = (1,1,1,0), so either (Qi,Qj) G" Aq 
or (Qk, QI) ^ Aq - a contradiction. 

Proof of claim (b) For a set S C Q denote [S] = ([i G S], [j G S], [k G S], [I G S], [k' G 
S], [I' G S]). Arcs (Qi,QV) and (Qk,QV) are not in Aq before sending the flow, therefore there 
exist sets S and S' with [S] = (1, ?,?,?,?, 0), [S] = (?, ?, 1, ?, ?, 0) and f Q (S) = f Q (S') = 0. We 
have (Qi,Qj),(Qk,Ql),(Qk',Ql') G Aq, therefore [S] = (1, 1, ?, ?, 0, 0), [S'} = (?, ?, 1, 1, 0, 0). 

Consider set S" = S U S" with [S"] = (1,1,1,1,0,0). Sets S and S' are minimizers of a 
submodular function /, and thus so is S". We have f Q (S") = f Q (S") - [is S"] + [j £ S"} = 
— 1 + 1 = 0, which implies the claim. 

Appendix B: proof of proposition [6] 

Let ip° be the input flow to the A-phase, S to be the set of nodes in V reachable from s in (N, A 2 ^) 
and (p = AdjustFlow((^°). Let c s i,cu and Jq be residual capacities and functions with respect to 
flow ip, and c° i5 c? t be residual capacities with respect to flow tp° . When the previous 2A-phase 
terminated, there were no augmenting paths from s to t in (N,A 2A ), hence A 2A cannot have arcs 
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(t, t) for i £ S and (s, i) for i G V - 5. Therefore, c° ti < [2A] — 1 for i £ 5 and c°, < [2A] - 1 for 
i E V - S. Define 

f A (S) = valued) +^>t+ J] 5 S i+^/^(5nQ) (20) 

iG5 iey-s QgQ 

Each augmentation in the A-phase preserves this equality (assuming that c s i , cu and Jq are up- 
dated accordingly). All residual values stay non- negative, therefore the number of augmentations 
cannot exceed (f A {S) — value{ip)) / [A]. Using ([20]) and the definition of step SO, we can write 



f A (S)- valued) < ]Tc° t + £ c° sl + £ 



F Q(5nQ) + ^|v9Q 4 -^| 



< n.([2Al-l) + ^a Q .rA]< 2n + £ a Q ) • [A] 

QgQ \ QgQ 

Appendix C: Proof of proposition [7] 

Proof of part (a) Let be distinct nodes in Q and (Qi,Qi'),(Qi',Qi") € ^4q. If 

Z{ < Zi" then obviously (Qi,Qi") 6 j4q. Suppose z» > in order to show (Qi,Qi") € j4q, 
we need to prove that mim-g^)^^//)^] p(/c) > 3A/2. Value Zj/ falls in one of the three intervals 
[.Zi,+oo), (— oo,Zj"], (zj»,Zj). These three cases are considered below. 

• Zi' > Zi > Zi". Since arc (Qi',Qi") belongs Aq and > Zi", we must have 

min g(k) > 3A/2 

ke[L(i')..R(i")-l] 

The claim then follows from the fact that L(i) > L(i') and so [L(i)..R(i") - 1] C [L(i!)..R(i") - I). 

• z% > Zi" > Zi'. Since arc (Qi,Qi') belongs to Aq and Zi > we must have 

min g(k) > 3A/2 
ke[L(i)..R(i')-l] 

The claim then follows from the fact that R(i') > R{i") and so \L{i)..R(i") - 1] C \L(i)..R{i') - 1] 

• Zi > Zi> > Zi". We must have 

min g(k) > 3A/2 min g(k) > 3A/2 

k£[L(i)..R(i')-l] ke[L(i')..R(i")-l] 

The claim the follows from the fact that R(if) > L(if) so [L(i)..R(i') - 1]U [L(i')..R(i") - 1] = 
\L(i)..R(i") — 1]. 

Proof of part (b) The transitivity of j4q and minimality of P implies that if (Qi, Qj) G P then 
the previous and the next arcs of P are respectively (i,Qi) and (Qj,j). The triple of consecutive 
arcs (i,Qi), (Qi,Qj), (Qj,j) will be denoted as and we will refer to it also as an "arc". Let 

Pq = (h, ji), ■ ■ ■ , {id, 3d) be the sequence of all such arcs of P (given in the order that they appear 
in P). Due to the minimality of P all 2d nodes involved must be distinct. It suffices to prove the 
proposition in the case when z% > Zj for all arcs {i,j) in this sequence. Indeed, if there are arcs 
{i, j) with Zi < Zj then we can push flow through them afterwards - as discussed in section 15.21 
this cannot violate the base polyhedron constraint. 
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We thus assume that Z{ > Zj for arcs £ Pq. Let and be two consecutive arcs 

in the sequence. We claim that Zj > z%i. Indeed, since path P is minimal, arc (Qi,Qj') is not in 
Aq. If Zii > Zj then (Qj,Qi') £ -Aq, so by transitivity we have (Qi,Qf) £ Ag - contradiction. 
If = then (Qi,Qi') £ Ag (since (Qi,Qj) £ Ag and z^ = Zj), so by transitivity we have 
(Qi, Qf) £ Aq - contradiction. 

We showed that z% t > Zj 1 > . . . > zi d > zj d . This implies that L(i\) < R(ji) < . . . < L(id) < 
R(jd)- Now consider k £ [0..m]; we need to show that g(k) = g{k) — J2k'=i z k > where z is the 
vector after sending [A] units of flow through P and g(-) is the corresponding residual function. 

It follows from the definition of z k that Q can be partitioned into two disjoint subsets S, T with 
k and m — k nodes, respectively, such that z% > z k > Zj for any i £ S, j £ T. Let us introduce the 
following terminology. Arc in Pq will be called left-exterior if Zj > Z, > z k + [A] (and thus 
i,j £ 5), and right- exterior if z fc — [A] > z.; > Zj (and thus i,j £ T). Clearly, after the update 
we have z% > Zj > z k for left-exterior arcs and z fc > Zj > Zj for right-exterior arcs (i,j). 

An arc in Pq is called exterior if it is either left-exterior or right-exterior, and interior otherwise. 
Note that an interior arc (i,j) must satisfy z% > z k > Zj, which is equivalent to the condition 
k £ [L(i)..R(j)]. This implies that Pq can have at most one interior arc. 

We now consider three possible cases. 

• All arcs in Pq are exterior. Then after the update we have Zj > z k > Zj for any i £ S, j £ T, 
so S contains k nodes % with the largest values of Zj. This implies that g(k) = g(k) — z(S). 
Since each arc in Pq either has both endpoints in S or both endpoints in T, we have 
z(S) = z(S), so g(k) = g{k) > 0. 

• Pq has an interior arc (u,v) with k £ [L(u)..R(v) — 1]; thus, <?(&) > 3A/2 since (Qu,Qv) £ 
Aq. We can assume without loss of generality that u £ S and « £ T. (Sets 5 and T could 
have been chosen in this way since L(u) < k and i?(v) > A;). After the update we have 
%i > -z 4 "' > Zj for any i £ <S, j £ T, so 5 contains fc nodes i with the largest values of Zj. This 
implies that g(k) = g(k) — z(S). Arc (u, v) is the only one in the sequence which has exactly 
one endpoint (namely u) in S. Therefore, z(S) = z(S) + [A] (where "+[A]" term comes 
from the update z u = z u + [A]), so g(k) = g(k) - [A] > 3A/2 - [A] > 0. 

• Pq has an interior arc (u, v) with R(v) = k. We must have u,v £ S and z u > z v = z k . After 
the update we have z v = z k — [A] and z» > z k > 2j for any i £ S — {v}, j £ T. Let (u',v') 
be the arc in Pq that immediately follows (u,v); if (u, v) is the last arc in Pq then we say 
that (u',v') does not exist. Two cases are possible: 

— Arc (u',v') does not exist or z k — 2 [A] > z u >. Then z v = z k — [A] > Zj for any 
j £ T. Thus, S contains k nodes i with the largest values of Zj. This implies that 
g(k) = g(k) — z(S). Since each arc in Pq either has both endpoints in S or both 
endpoints in T, we have z(S) = z(S), so g{k) = g{k) > 0. 

- Arc (vf, v') exists and z u < = z k - [A] ; thus, L(v!) = R(v) + 1 = k + 1, z k+1 = z k - [A] . 
After the update z v = z k — [A], z u > = z k , so the set S' = (S — {v}) U {u'} contains k 
nodes i with the largest values of Z{. This implies that g(k) = g(k) — z(S'). We have 
z(S') = z(S) + [z u i — z v ] = z(S) + [A], so g(k) = g(k) — [A]. We now need to show 
that g(k) > [A]. 

Conditions (Qu,Qv), (Qu' ,Qv') £ A^ imply that g(k - 1) > [3A/2] and g(k + 1) > 
[3 A/2] . We can write 

g (k)-g(k-l) = [g(k) - g(k - 1)] - z k 
g(k + l)-g(k) = { g (k + l)-g(k)]-z k + 1 
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Since g(-) is concave, we have g(k) — g(k — 1) > g(k + 1) — g(k). Thus, 
g(k) - g(k -l) + z k > g(k + 1) - g(k) + z k+1 



2g(k) > g(k-l)+g(k + l)-[z k -z k+1 ] 

> [3 A/2] + [3 A/2] - [A] > 3 [A] - [A] = 2 [A] 

Proof of part (c) Let S be the set used in eq. (fl"4"|) . and denote T = Q — S, k = \S\. We assume 
that S ^ and S ^ Q, otherwise the LHS in eq. (fl"4"|) would be 0. Let i be a node in S with 
the minimum value of zi and j be a node in T with the maximum value of Zj. Since there was no 
augmenting path upon termination of the previous 2A-phase, set cannot have arc (Qi,Qj)- 
Therefore, Z{ > Zj and g(k) < [3AJ for some k € [L(i)..R(j) — 1]. The choice of i,j and condition 
Zi > Zj imply that min{zj/ j zi> € S} > max{zj/ | zy € T}, hence fQ(S) = g(\S\) = g(k). Thus, 
we need to show that g(k) < olq ■ [A] where cxq = 3(m — 1). If k = k then the claim is obvious. 
Suppose that k ^ k. Two cases are possible: 

• k > k. We have k + l<k + l< R(j), so z k+1 > z R ^ = zj. We cannot have z k+1 > zj 
since in this case there would be at least k + 1 nodes f £ Q with Zj> > zj; by the choice of 
j these nodes would belong to 5, so we would have \S\ > k + 1 - contradiction. Thus, we 
must have z k+1 = z k+2 = . . . = z R ^\ This implies that function p(k') = Y^,k"=i zk " ls nnear 
in [k..R{j)\. We have g(k') = g(k') — p(k') where g(-) is a concave function, therefore g(-) is 
also concave in [k..R(j)]. There holds k € [k + l..R(j) — 1], thus 

We have g(k) < [3A\ and g(R{j)) > 0, therefore g(k) < R{j) ~ k [3AJ < (m - 1) • L3AJ < 

R(j)-k 

• k < k. We have L(i) < < A;, so Zj = z L ^ > z k . We cannot have Z{ > z fc since in this 
case there would be at least m — k + 1 nodes i' € Q with ^ > Zi> ; by the choice of i these 
nodes would belong to T, so we would have \T\ > m — k + 1 - contradiction. Thus, we must 
have z L ® = . . . = z k ~ l = z k . This implies that function p(k') = X^"=i zk " ^ s linear in 
[L(i) — l..k]. We have g(k') = g(k') —p(k') where g(-) is a concave function, therefore g(-) is 
also concave in [L(i) — l..k]. There holds k € [L(i)..k — 1], g(k) < [3A\ and g(L(i) — 1) > 0, 
so similar to the previous case we conclude that g(k) < fc ~ (L(i) ~ 1} [3AJ < (m - 1) • [3AJ < 

<xq- [A]. 



Appendix D: Implementation of GetNeighborsg(Qz) for cardinality- 
dependent terms 

Case 1 Assume that fo(S) = g(\S\) for S Q Q where g is concave. We use the same notation 
as in section [531 

First, let us describe the data structure for Q. Nodes i € Q will be grouped into "supernodes" 
according to their value of The set of supernodes is denoted as Q. The cardinality of Q equals 
the number of unique values in the set {zi \ i S Q}. At each supernode u € Q we store values 
z u = Zi, L(u) = L(i), R(u) = R(i) where i is a node contained in u. We treat supernode u as 
the set u = {Qi | i 6 Q,Zi = z u }. Supernodes u sorted by their values of z u will be stored in a 
doubly-linked list. Each u S Q also have a pointer to a doubly-linked list of nodes in u, and each 
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node Qi G N will have a pointer to u G Q with Qi G u. Finally, we maintain residual function g(-) 
as an array of size 0(m). It is easy to see that after each augmentation this data structure can be 
dynamically updated in 0{m) time. 

For each supernode u we maintain flag REACHED(u) = /\ i£u REACHED(<5i); at the beginning of 
the BFS it is set to false. Procedure GetNeighborSg((5i) is defined as follows: 
GetNeighbors^(Qi) 

• Set REACHED(Qi) := true and S := 0. Determine supernode u with Qi G u. 

• If REACHED(u) is true then stop, otherwise set REACHED(u) := true and continue. 

• If min fc6 [ L ( u )..fl( u )_ 1 ] g(k) > 3A/2 call Add(u). 

• If u has left neighbor u~ with z u - > z u call Add(u~) and ProcessLef t(u~). 

• If u has right neighbor u + with z u > z u + and KU&ke[L(u)..R(u+)-l] 9~(k) — 3A/2 call Add(n + ) 
and ProcessRight(u + ). 

ProcessLef t(u) 

• If REACHED(n) is true then stop, otherwise set REACHED(ii) := true and continue. 

• If u has left neighbor u~ with z u - > z u call Add(u~) and ProcessLef t(u~). 
ProcessRight(n) 

• If REACHED(u) is true then stop, otherwise set REACHED(ii) := true and continue. 

• If u has right neighbor u + with z u > z u + and ^^ke[L<u)..R(u+)—l] 9~{k) — 3A/2 call Add(n + ) 
and ProcessRight(u + ). 

Add(u) 

• For each node Qi G u with REACHED(Qi) = false set REACHED(Qi) := true and add Qi to 
S. 

The correctness of this procedure should be clear. Note, ProcessLef t(u) and 
ProcessRight(ti) are only called when some node Qi € u has been reached by BFS. If REACHED(u) 
is true then all nodes that can be reached from Qi (and from other nodes in u) via arcs in Aq have 
already been added, which justifies statement "If REACHED(u) is true then stop". Steps following 
this statement will be executed at most once for each supernode u, therefore each node, supernode 
and element of array g(-) is accessed at most constant number of times during a single BFS search. 
Thus, Pq = 0(m). 

Case 2 We now assume that /q(5) = g(\S n Q'\, \S D Q"\) for S C Q where Q' , Q" are disjoint 
subsets of Q. Without loss of generality we can assume that Q = Q' U Q". Denote ml = \Q'\, 
m" = \Q"\, m = \Q\ = ml + m" . Let y € and z € be vectors with yi = ipQi for i G Q' 
and Zj = ifQj for i G Q" . We define sequences (y , ...,y m ) and {z l ,... 7 z m ) similar to the 
case above; y k and z k are the fc-th largest numbers among values y% and Zi, respectively. Indexes 
L(i) and R(i) are also defined as in section I5.2( we have 1 < L(i) < R(i) < ml for i G Q' and 
1 < L(j) < R(j) < ml 1 for j G Q" . Let g(k', k") = min{/ Q (5) | \S n Q'\ = k' , \S n Q"\ = k"} be 
the "residual function". We have 

k' k" 
g{k\k") = g{k' 1 k")-Y j y k -Y, zk 

k=l k=l 
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It can be seen that <?(■,•) is a Monge matrix [3J, i.e. for any < k[ < k! 2 < m' and < k'{ < k 2 < m" 
there holds g(k[,k'{)+g(k' 2 , k%) < g(k[, k%)+g(tf 2 , k'{). This follows from f Q (S'nS")+f Q (S'uS") < 
/q(S') + fc}{S") where 5" contains first nodes of Q' and first ^ nodes of Q" , and 5" contains 
first k 2 nodes of Q' and first k'[ nodes of Q". (We assume that nodes in Q' and Q" are sorted so that 
components yi and Z{ are non- increasing.) For a row fc' G [0..m'] let k"(k') G [0..m"] be the column 
that contains the leftmost minimum entry in row k' . Thus, k"(k') = minjfc" G [0..m"] | g(k' , k") = 
miiifc// g(k', k")}. It is known [1] that Monge matrices are monotone, i.e. k"(0) < k"(l) < . . . k"(m'). 
Furthermore, they are totally monotone, i.e. every submatrix is monotone. As shown by Aggarwal 
et al. pQ, indexes k"(0), . . . , k"{m') for a totally monotone matrix can be computed in 0{m) time. 

We can describe data structures for implementing GetNeighborsg(Qi). Nodes in i G Q' will be 
grouped into supernodes according to the values yi analogously to case 1. A similar data structure 
will be used for nodes in Q" . We will maintain an array of cumulative sums X^fc=i V k f° r ^' ^ [0..m'] 
and X]fc=i zk f° r k" G [0..m"], which will allow computing g(k' , k") in O(l) time. At the beginning 
of each BFS we will compute indexes k"(k') for k' G [0..m'] using the algorithm in [I] and also 
indexes k'(k") = min{A/ G [0..m'] | g(k' ', k") = min^/ g(k' , k")} for each column k" G [0..m"]. 

Arcs in Aq can be split into four groups ^4ocb -^lO) ^11 where A a p = {(Qi,Qj) G Aq \ [i G 
Q'\ = a, [j G Q"\ = /?} and [•] is the Iverson bracket: it is 1 if its argument is true, and otherwise. 
Consider the version of GetNeighborsg(<5i) that processes only arcs in a specific group, rather 
than all arcs in Aq. It suffices to show how to implement such procedure for each of the four 
groups; these procedures will be called sequentially. 

First, consider arcs in ^n. Using the same argumentation as in section l5\2l we conclude that 
sending flow from a node Qi to another node Qj (i,j G Q') is possible if and only if one of the 
two conditions hold: (a) yi < yf, (b) y, t > yj and mmk'e[L(i)..R{j)-i]9'(k') > 1 where we defined 
g'(k') = minfc// g(k', k"). Thus, the set ^oo is constructed completely analogously to the set Aq 
in section T5.2I except that function g is replaced with g' and threshold 3A/2 is replaced with 1. 
Accordingly, we can use an obvious adaptation of the procedure for the case 1. Note, g'(k') can 
be evaluated in O(l) time using arrays of indexes k"(k') and cumulative sums for vectors y and z. 
Arcs in ^oo can be handled in a similar way. It remains to show how to handle arcs in A±q (the 
set Aqi will follow by symmetry). 

Consider nodes i G Q' , j G Q". Sending [A] units of flow from i to j via Q, i.e. the operation 
Hi '■= Hi + [A], Zj := Zj — [A], affects function <?(•,•) as follows: values g(k',k") for (k',k") G 
[L(i)..m'] x [l..R(j) - 1] are decreased by [A] and values g(k',k") for (k',k") G [l..L(i) — 1] x 
[R(j)..m"] are increased by [A]. Thus, sending flow is possible if and only if 



mm 



{g(k', k") | (k' , k") G K(L(i), R(j))} > , K(a, b) = [a..m ; ] x [0..6 - 1] 



There holds K(a, b\) C K{a, 62) for 61 < hi, therefore the set of arcs in A\q from Qi have the form 
{(Qi,Qj)\R(j)<b(L(i))} where 



b(a) = max < b G [l..m"l min q(k' , k") > > a G [l..m'l 

I (k',k")eK(a,b) V J 



(21) 



(If the set in (j2T|) is empty then we assume that 6(a) = 0.) We compute indexes 6(a) at the 
beginning of BFS in linear time using the following recursion: 

b(m ) = k (m ) (since min^" g(m', k") = g(m',m") = 0) 

6(a) = ( 6(a + 1) if^»)>0 W ae[l..m>-1] 

\min{6(a + l),A;"(a)} if g(a, k"(a)) = 

Note that < 6(1) < ... < b(m') < m" . Procedure GetNeighbors^(Qz), % G Q' for the set of 
arcs Aiq is implemented as follows. First, we locate the rightmost supernode v C Q" satisfying 
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R(v) < b(L(i)). (Pointers to these supernodes for each supernode u C Q' can be computed at the 
beginning of BFS.) We then call procedure Add(u), which is defined as in the case 1, and procedure 
ProcessLef tio(f ) defined as follows: 
ProcessLef tio(v) 

• If REACHED(v) is true then stop, otherwise set REACHED(t> ) := true and continue. 

• If v has left neighbor v~ with z v - > z v call Add(t;~) and ProcessLef t\o(v~). 
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